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SUMMARY 

. Vibration analysis of pretwisted blades of asymmetric airfoil cross sec- 
tion is performed by using two mixed variational approaches, one proposed by 
Dean and Plass, and the other by Reissner. Numerical results obtained from 
these two methods are compared to those obtained from an improved finite- 
difference method developed by the authors, and also to those given by the 
ordinary finite-difference method. The relative merits, convergence properties 
and accuracies of all four methods are studied and discussed. The effects of 
asymmetry and pretwist on natural frequencies and mode shapes are investigated. 
The improved finite-difference method is shown to be far superior to the con- 
ventional finite-difference method in several respects. Close lower bound 
solutions are provided by the improved finite-difference method for untwisted 
blades with a relatively coarse mesh while the mixed methods have not indicated 
any specific bound. 


INTRODUCTION 

Turbine and compressor blades, helicopter rotor blades and propeller 
blades are generally pretwisted and possess asymmetric airfoil cross sections. 
For such blades, the centroid, center of flexure and center of torsion are 
noncoincident. General displacements of an asymmetric blade, even without 
pretwist, consist of translations coupled with rotations. Duncan, Ellis and 
Scruton (ref. 1) showed that for long members having rigid or quasi-rigid sup- 
port, the center of flexure and torsion center are very nearly coincident. 
Carnegie (ref. 2) extended their work for asymmetric airfoil blades and deter- 
mined the coordinates of the center of flexure with respect to the centroid of 
the blade cross section. When a blade of asymmetric airfoil cross section 
vibrates, coupled bending-torsion vibrations occur. Coupling between the bend- 
ing motions in the two principal directions will be intensified due to pretwist 
of the blade. Further, an increase in torsional rigidity, over and above that 
due to St. Venant, takes place due to pretwist. For blades having low aspect 
ratios and a wide range of thickness ratios, the effect of warping rigidity is 
very important (ref. 3). 

The foundations of the coupled flexure-torsion theories are based on the 
thin-walled open section formulation and date back to 1936 when Wagner (ref. 4) 
derived the relation between twisting moment and torsional rotation for the 
case of nonuniform torsion. The coupled bending-torsion equations of motion 
were derived by Garland (ref. 5), Gere and Lin (ref. 6), Houbolt and 
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Brooks (ref. 7), Montoya (ref. 8) and Yu (ref. 9), among others, using the 
classical Euler-Bernoul 1 i theory allowing for various complexities. Using the 
Timoshenko beam theory, which allows for shear deflection and rotary inertia 
effects, and incorporating further complexities, the coupled equations of 
motion were derived by Carnegie (ref. 10), Rao and Rao (ref. 11) and 
Subrahmanyam, Kulkarni and Rao (ref. 12) among others. Recent research in hel- 
icopter and turbo-prop blades indicates that the linear equations of motion are 
inadequate for establishing flutter boundaries. Geometric nonlinear theory is 
to be used with various degrees of complexity. Some such theories were pre- 
sented by Hodges and Dowell (ref. 13), Rosen and Friedmann (ref. 14), and Kaza 
and Kvaternik (ref. 15) among others. 

Several methods of solution of vibration problems are available. Broadly, 
these are classified as belonging to either the continuum model approach or to 
the discrete model approach. A comprehensive review of the work done using the 
various methods in blade vibration area can be found in the review articles by 
Rao (ref. 16) and Leissa (ref. 17) among others. 

An examination of the published literature reveals that among the various 
methods used for solution of the differential equations of motion, the first- 
order finite-difference method has attracted the greatest attention while 
relatively few works exist which use higher-order finite-differences. Even in 
those few works using second-order central differences either the approxi- 
mations made to eliminate fictitious stations were not good enough to give 
accuracy superior to first-order theory results (ref. 18) or a method of using 
forward or backward differences at beam boundaries was used to avoid the ficti- 
tious stations (ref. 19). An improvement to the second order finite-difference 
method was made by the present authors and this refined method was applied to 
uncoupled vibration analysis of tapered cantilever beams in (ref. 20), and to 
pretwisted tapered blades executing coupled bending-bending vibrations in 
(ref. 21). However, the applicability of the method for vibration analysis of 
pretwisted asymmetric blades has not been investigated. 

Another area, currently receiving considerable attention, is the utilization 
of mixed or hybrid variational techniques for vibration analysis. Among these 
mixed variational principles, the Reissner method (ref. 22) has received the 
greatest attention. For solution of the dynamic vibratory problems, the usual 
way of formulating the dynamic principle is to consider the difference of the 
kinetic energy and the Reissner functional, consisting of contributions from 
potential and complementary energies together with the work done by the body 
forces and surface tractions (ref. 23), and to minimize the same according to 
the Ritz process. One thus has the flexibility to vary displacements and 
stresses independently and can obtain simultaneously good distributions of all 
the independent parameters. Another more flexible way of formulating the 
dynamic variational functional is demonstrated by Dean and Plass (ref. 24) 
where the kinetic energy effects were incorporated by using linear momentum 
density in conjunction with the Reissner functional. This way to include 
angular momentum and rotary inertia effects was demonstrated by the present 
first author's earlier work (refs. 25 and 26), and it was observed that both 
mixed variational principles (refs. 22 to 24) lead to identical results if 
identical shape functions and beam kinematics were used in the analysis. How- 
ever, attempts to incorporate coupling between the flexural and torsional 
degrees of freedom were not made so far. An important advantage offered by 
the Dean and Plass method is that one can handle the gyroscopic forces in the 
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vibratory problem in a straightforward manner. Since this variational prin- 
ciple contains products of linear or angular momenta and the corresponding 
velocities in the kinetic energy equivalent, and since the gyroscopic forces 
also are functions of velocities, the variational functional possesses consist- 
ently first order time derivative functions. Unlike the Reissner or potential 
energy methods, which require either a transformation procedure or other more 
complicated methods for solution {refs. 27 and 28), the Dean and PI ass method 
leads directly to a standard eigenvalue problem. 

The objectives of the present work are: (1) to apply the first, and second 

order finite-difference methods, the Dean and Plass method and the Reissner 
method for the general case of vibration with coupling between bending and 
torsional degrees of freedom; (2) to assess the computational advantages and 
disadvantages of each method; and (3) to study the effects of pretwist and 
asymmetry in coupling the component modes. Numerical results from all four 
methods will be presented and comparisons will be made to those from other 
existing theoretical and experimental results (refs. 29 to 33). 
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SYMBOLS 

area at any section 
matrices 

arbitrary parameter in shape functions 

body force distribution 
torsional stiffness 
warping rigidity 

center-of-f lexure of blade cross section 
Young's modulus 
shape functions 
modulus of rigidity 
length of each elemental beam segment 
second moment of area of a cross section about x-^ 
axis 

second moment of area of a cross section about y^y^ 
axi s 

product of inertia of a cross section about x^x^ and 
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Dean and Plass dynamic variational functional 
dynamic Reissner functional 

time averaged values of Ip and Ir respectively 
dummy indices 

number of terms in the assumed solutions 

length of blade 

bending moments 

number of beam segments 

natural radian frequency 

coordinates of center-of-f lexure with respect to centroid at 
any section distant n from root 

r x /L, r y /L etc. 

coordinates of center of flexure with respect to centroid at 

the root section . : 

polar radius of gyration, y{(I xlxl + X y i y i) /A + r x + r y } 
part of boundary over which surface fractions are prescribed 

surface fractions 
time 

kinetic energy 
twisting moment 
displacement field 

displacement in x, y, and z directions 
volume 

components of linear momentum per unit volume 
linear momentum density along x-^x^ and y^ 
directions, respectively 
angular momentum density 

displacement of center-of-f lexure in xx, yy directions 
respectively 

displacements of centroid in x^x^, y^y-^ directions 
respectively 

coordinates measured with respect to center of flexure and 
centroid respectively 
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coordinates axes through centroid and center of flexure 
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longitudinal axis 

coordinate distance measured along the length of blade from 
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mass density 
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DEVELOPMENT OF FREQUENCY EQUATIONS 

Frequency Equations by Finite-Difference Methods 

For a pretwisted blade of asymmetric airfoil cross section having a total 
pretwist angle y over length L, the equations of motion (A18), appendix A, 
assume the following nondimensional form: 
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The coefficient a^, b^ ... h]_, a 2 , b 2 ... h 2 , a 3 ... d 3 are all functions 
of n and are presented in appendix B. 
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In implementing a finite difference procedure for the solution of the 
equations of motion, one substitutes the finite-difference expressions for the 
derivatives in the differential equations and eliminates the fictitious sta- 
tions out side the beam domain by enforcing the boundary conditions. The 
resulting equations contain the displacements x-j, y-j and ei at an arbitrary 
station i. The coupled equations are then evaluated at each station of the 
beam divided into n segments, i = l,2...n. One thus obtains a set of 3n 
simultaneous equations. 


In references 20 and 21, the present authors developed a refined procedure 
for eliminating fictitious stations that arise in using second order central 
differences. Use is made of the recursive relations and the finite-difference 
expressions reported in reference 20 in the present investigation, and the 
resulting equations are represented in the familiar form of the eigenvalue 
problem 
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In the preceding equation, &, £, £,, Q and £ are all square matrices of order 
(nxn), {x-j}, {y-j} and { } are column matrices containing the linear and 
torsional displacements of the n-stations and F, G, H, I, J, K and L are 
diagonal. The matrix Q is a null matrix. Further, each submatrix A, B, £ 
or D has a band width of five for the first-order theory and seven for second- 
order theory. The submatrix £ has a band width of three for first-order theory 
and five for second-order theory. The matrices fl, B, £, and £ are 
nonsymmetric in first and second order finite-difference applications. For 
brevity elements of these matrices are not presented, but one can easily develop 
these matrices with the help of references 20 and 21. 


Frequency Equation from the Dean and Plass Principle 

The Dean and Plass dynamic variational functional is developed in 
appendix A and has the functional relationship of the form 

Ip = ^D ( x »y » P x 1 ’ Pyl * P 9* ^x 5 ^y * "^" 9 * ^ ^ ) (5) 

The parameters in the functional given by equation (5) are independent of each 
other and depend only on n and t. 

The following shape functions in series form (refs. 23 and 25) are assumed 
for solution of the coupled bending-bending-torsion vibration problem: 
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The parameters Rj, R^ . . . R5, G(t), H(t) etc. are defined in appendix B. 
The shape functions given in equation (2) satisfy the boundary conditions 
(ref. 12) applicable to a cantilever beam fixed at n = 0 and free at n = 1. 
Substituting equations (6) and their appropriate derivatives into equation 
(A13), eliminating the time dependence according to the averaging procedure 
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and applying the Ritz process to minimize Ip with respect to the arbitrary 
parameters, that is. 
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one obtains the frequency equation which can be written as 

A + P n B = 0 
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In equation (9), & and B are symmetric square matrices. For brevity elements 
of these matrices are not presented here. 
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Frequency Equation from the Reissner Method 
The Reissner variational functional (ref. 22) is given by 
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Neglecting body forces, surface tractions, shear deflection, rotary inertia, 
thermal and rotational effects, one can show that (ref. 12) 
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Using similar shape functions as were used for the Dean and Plass method and 
proceeding on the same lines as described earlier, one can obtain the frequency 
equation which can be written as 


C + P 2 D = 0 
n~ 


( 12 ) 


Matrices C and £ are symmetric, the elements of which are not presented here 
to save space. 


METHOD OF SOLUTION 

The eigenvalue problems given by equations (4), (9) and (12) were solved 
by using a standard eigenvalue extraction routine. This library routine deter- 
mines all the eigenvalues and the associated eigenvectors. For the eigenvalue 
problem formulation using mixed variational approaches, certain integrations 
were to be performed numerically. This was accomplished by using a 15 point 
Gaussian quadrature formula. The computer programs developed in FORTRAN lan- 
guage were run on an IBM 370 computer at the NASA Lewis Research Center. The 
lowest ten coupled frequencies and associated mode shapes were determined for 
studying the effects of asymmetry and pretwist in coupling the modes using 
typical airfoil blade data (ref. 2). Results are presented and discussed in 
what follows. 


RESULTS AND DISCUSSION 

In order to study the accuracies and computational efficiencies of the 
four methods discussed in this report, airfoil blade data developed by Carnegie 
(ref. 2) and used by Carnegie and others for numerical calculations (refs. 29 
to 32), was adopted. Such numerical data pertaining to the typical blade pro- 
file shown in figure 1 is presented in tables I and II. It may be noted that 
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the coordinates of the center-of-f lexure measured with respect to the centroid 
assumed in this report are the experimental values (ref. 2). Further, the 
values of torsional rigidity used in the present calculations are taken from 
the experimental data presented by Carnegie (ref. 29). By using the numerical 
data presented in tables I and II, results obtained from the four methods are 
compared mutually and to those existing in the literature in what follows. 


Relative Convergence Rates and Accuracies 

A typical set of results obtained by using first and second-order central 
difference formulations and with various mesh sizes are presented in table III. 
Results given by both first and second order theories are lower bounds for an 
untwisted blade. No specific bound is given by the finite-difference methods 
for twisted blade cases (ref. 21). 

From the results presented in table III, one can state that the conver- 
gence shown by the second order theory is much more rapid than the first. 
Further, the accuracy of the converged results can be assessed by comparing 
the natural frequencies given by the finite-difference solution with a given 
mesh size (number of beam segments) to corresponding extrapolated results using 
the Richardson extrapolation procedure (ref. 34). For this purpose, the 
results using n = 10, 20 and 30 in the second order theory are extrapolated 
and are presented in the last column of table III. In comparison with the 
extrapolated result, the first-order theory with n = 60 shows maximum errors 
of the magnitudes of about 0.4 percent and 0.64 percent in the lowest seven 
coupled mode frequencies for the cases y = 0° (untwisted blade) and y = -30° 
(clockwise pretwisted blade) respectively. Corresponding maximum errors from 
the second order theory are 0.075 percent and 0.322 percent when the beam is 
divided into 30 segments (n = 30). 

From these observations, it is concluded that the second-order theory 
produces results with accuracies of practical interest with a coarse mesh size. 
Results given by the second order theory with n = 30 are superior to those 
from the conventional first order theory with n = 60 and thus, the extrapola- 
tion procedure is not necessary when the improved second order theory is used 
with a suitable mesh size. 

In order to assess the relative merits of the two mixed variational prin- 
ciples namely the Reissner method and the Dean and Plass method, identical 
shape functions were assumed for the independent variables in both methods. 
Natural frequencies were determined for the two typical values of pretwist 
(y = 0° and y = -30°) considered earlier while using finite-difference 
approaches. Natural frequencies were obtained by varying the number of terms 
in the assumed solution for each shape function. The convergence shown by the 
methods was found to be almost identical. These results are shown in 
tables IV(a) and (b). This is an expected trend since the momentum functions 
were chosen to have identical shape functions as those corresponding to the 
respective displacement functions, thus, making the two variational methods 
numerically equivalent. However, the flexibility of independent variations of 
displacements and momenta can effectively be utilized by incorporating differ- 
ent forms of shape functions as shown in reference 24 and 26, consistent with 
the boundary conditions. 
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From the identical convergence pattern produced by both mixed variational 
principles, the accuracy of the theoretical development of the Dean and Plass 
functional is conclusively verified. A further comparison of the results pre- 
sented in table IV to the corresponding sets presented in table III indicates 
that the lowest seven coupled modes obtained by using a six-term solution of 
the Dean and Plass or Reissner method agree very closely with the extrapolated 
results from finite-difference solutions. Furthermore, the mixed variational 
principles considered here have not shown any specific bound, and the results 
presented in table IV can be seen to have an oscillatory nature. 

Finally, a comparison of the present converged results is made to the 
theoretical and experimental results available in the literature. Such a com- 
parison is o shown in figure 2. Both clockwise and anticlockwise pretwists rang- 
ing from 0° to 89° were considered. Results produced by finite-difference 
procedure using second-order central differences with n = 30 and those pro- 
duced by the Reissner and, Dean and Plass methods with K = 6 were close for 
the lowest seven modes. These results agree very closely with the theoretical 
results of Carnegie et al. (ref. 30) obtained by using a transformation tech- 
nique and also with experimental results. Although not shown here, the mode 
shapes obtained by using the various methods were also close. 


Effects of Pretwist and Asymmetry in Coupling the Modes 

For blades of doubly symmetric cross section, the effect of pretwist in 
coupling the two principal bending motions has been well understood in the 
published literature. Rosard (ref. 35) observed the effect of pretwist on the 
coupled bending-bending frequencies for the first time. Later, Carnegie et al. 
(refs. 29 to 31), Rao (ref. 16) and the present authors (refs. 3, 21, and 23), 
among others, studied the effect of pretwist on the natural frequencies and 
mode shapes. However, the combined effects of asymmetry and pretwist on the 
natural frequencies and mode shapes have not been completely understood. 

In order to understand any specific coupling trends that may exist in the 
coupled bending-bending-torsional frequencies of pretwisted blades of asym- 
metric airfoil cross section, it is necessary to establish the coupling trends 
for simple cases of blade geometry first. For this purpose, using the second 
order finite-difference method with 30 beam segments, the uncoupled bending 
and torsional frequencies and mode shapes were determined by assuming r x0 = 0, 
r„ 0 = 0 and y = 0 in the general computer program. These natural frequen- 
cies are presented in table V under column 2. In this table, the frequencies 
of the beam vibration in the flexible plane are designated by FI, F2, ..., F5 
for the lowest five modes, those in the stiffer plane are designated by SI, S2 
for the lowest two modes while Tl, T2 and T3 represent the lowest three tor- 
sional frequencies in the increasing order. The coupled bending-bending fre- 
quencies and the associated mode shapes were next determined by successively 
assigning different pretwist angles. The coupled bending-bending frequencies 
for typical pretwist angles are presented in table V under columns 3 to 6. 
Although the mode shapes were determined for all the pretwists considered here, 
only one such set (pretwist angle of 89°) is presented in figure. 3. Using 
these uncoupled mode shapes or the coupled bending-bending mode shapes, one 
can classify any coupled bending-torsion or coupled bending-bending-torsion 
mode shape as belonging closely to a basic mode category. Further results 
were generated for various combinations of the coordinates of asymmetry and 
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pretwist angles. The modes were then classified as belonging to a basic mode 
category by examining the associated mode shapes in relation to either the 
uncoupled mode shapes or to the coupled bending-bending mode shapes. These 
results are tabulated in tables VI to VIII. Thus, in table VI, for blades with 
r x0 = 0.193 mm, r y0 = 0 are presented, in table VII those for r x0 = 0, 
r = 1.1938 mm are shown while in table VIII, the results of general asym- 
metric airfoil blade with r x0 = 0.193 mm, r y0 = 1.1938 mm are listed for 
various angles for pretwist. It may be noteo that the frequencies are not 


necessarily in an increasing order as the mode number is increased since the 
modes are classified as belonging to a basic mode category. Only one set of 
coupled bending-bending-torsion mode shapes are presented in figure 4 for the 
case y = 89° for brevity. The general coupling trends as observed from these 
results are presented below. 


(1) For pretwisted blades of doubly symmetric cross section, strong 
coupling exists between the second uncoupled flapwise mode (F2) and fundamental 
uncoupled chordwise mode (Si). As pretwist increases, the lower frequency 
decreases and the higher frequency increases to form the corresponding coupled 
frequencies. Similar coupling between the uncoupled modes S2 and F5 exists 
but the magnitude of frequency variations are more severe for the set (F2, SI) 
than for the set (S2, F5). This frequency variation trend is consistent with 
earlier observations reported in references 23 and 35. The coupled bending- 
bending mode shapes obtained in the present investigation are compared to those 
in reference 33. Consistent trends have been observed in all the results 
compared. 


(2) For untwisted blades with asymmetry about only one plane, coupling 
between bending in one principal direction and torsion occurs. Results for 
such a case with r y0 = 0 are shown in column 3 of table VI, and for the case 
r = 0 in table VIT. From the results presented in table VI, one can 
observe that the frequencies corresponding to (F2, Tl) and (F3, T2) couple in 
such a way that the lower uncoupled frequency in each pair is reduced and the 
higher one is increased to form the corresponding coupled bending-torsion fre- 
quencies. Although not shown here, similar coupling trends exist for the pairs 
(F4, T3) and (F5, T4). Similar coupling trends are observed for the coupled 
modes having coupling between chordwise bending and torsional motions. This 
coupling trend is shown by the frequency pair (SI, Tl) in table VII. Although 
not shown, the frequency pair (T3, S2) also shows similar coupling trends. 

(3) For pretwisted blades with asymmetry about only one plane, the cou- 
pling trends are combinations of those discussed in (1) and (2) above. A care- 
ful examination of the results presented in tables VI and VII for various 
values of pretwist angles in relation to the corresponding ones in table V 
indicates that a specific trend in the three close frequency values correspond- 
ing to F2, SI and Tl modes exists. The frequencies corresponding to the pre- 
twisted symmetric blade listed in table V will be further coupled because of 
asymmetry in such a way that the lower two frequencies will be reduced and the 
highest one will be increased to result in the coupled bending-bending-torsion 
frequencies. Other secondary coupling trends are observed for the frequency 
pair (T2, F4) when the frequencies are very close to one another, in table 
VII, and for all pretwist angles in table VI. The reason for this is that the 
flap-torsion coupling is more important for blades having a symmetry about the 
chord (table VI) than for those with asymmetry about the chord (table VII). 


11 



Similar coupling trends exist for the frequencies corresponding to the 
frequency triplet (T3, S2, F5). However, results for these modes are not pre- 
sented here because obtaining accuracy in these higher mode frequencies would 
have required a greater number of beam segments. 

(4) The coupling trends for pretwisted asymmetric airfoil blades present- 
ed in table VIII are similar to those observed for the simpler cases considered 
earlier. It has been observed that the coordinate ry 0 is more important . 
than r xo in coupling the component modes. 

(5) From the mode shapes obtained for all the cases studied in this 
investigation, it is observed that the coupled modes which are well separated 
from a basic torsional mode frequency are not seriously effected by the tor- 
sional coupling (refer figs. 3 and 4). These mode shapes can be closely 
approximated by the corresponding coupled bending-bending mode shapes. For 
modes which are closer to a basic torsional mode, the extent of torsional cou- 
pling must be determined. 

The trailing edge deflections calculated in this investigation show close 
agreement with those presented by Carnegie et al. (ref. 30). Although the 
trailing edge deflections were determined for all cases of pretwist angles 
studied, they are not presented in this report for brevity. 

Finally, the coupled frequencies were determined for different sets of 
asymmetry coordinates, (r xo = 2.54 mm, ry 0 = 1.1938 mm; and r x0 = 0.193 mm, 
ry 0 = 2.54 mm), to examine whether there are any changes in the coupling 
trends with large variations in these coordinates. Although these results are 
not presented here, it was observed that the general coupling trends remain 
unchanged. 


RELATIVE COMPUTATIONAL EFFICIENCIES 

In order to evaluate the relative computational efficiencies, the CPU 
time required by each method for a typical blade example was determined. These 
values are presented in table IX. 

Comparing the CPU times required by firstand second-order finite-difference 
methods for a given number of beam segments, one can state that both methods 
require nearly the same amount of CPU time. The maximum variation of CPU time 
between the two methods is of the order of +3 percent. If one considers the 
level of accuracy for estimation of computational efficiency, the first-order 
finite-difference method requires a solution with 60 beam segments to match the 
accuracy of results produced by the second-order finite-difference method with 
n = 30. The computation time required for first-order theory would be nearly 
13 times that required by the present improved second order theory. Finally, 
if one considers an extrapolation procedure as a suitable alternative, then, 
both first and second order theories presented here are capable of taking the 
advantage. Thus, it is preferable to use the present improved finite-difference 
procedure as compared to the classical first-order theory. 

When a comparison is made between the computational times required by the 
Reissner and the Dean and Plass methods, one observes that the Dean and Plass 
method requires nearly twice as much time as the Reissner method. For the 
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particular example considered here with identical shape functions in both 
methods, it has been observed that the convergence trends and accuracies given 
by these two methods are identical. Thus, for the class of vibratory problems 
treated here, it is advantageous to use the Reissner mixed variational princi- 
ples as compared to the Dean and Plass principle. On the other hand, for solu- 
tion of problems having gyroscopic forces (Coriolis accelerations), the usual 
formulation by the potential energy method or the Reissner method requires a 
transformation of the time dependent variables by a linear transformation to 
attain the eigenvalue problem in a standard form (ref. 28) while the Dean and 
Plass method seems to offer a more direct solution facility. This advantage 
together with any associated computational advantages relative to the existing 
approaches are yet to be assessed. 

When a choice is to be made between the finite-difference procedures and 
the energy methods for potential computational needs, it appears that energy 
methods have a great deal to offer when computational time is the principle 
constraint. As can be seen from tables III and IV, one can obtain the lowest 
seven coupled modes accurately with a six-term solution using the Reissner 
method. These results can be obtained in about 13 percent of the CPU time 
required by the second order finite-difference method. While making this com- 
parison, it may be noted that the symmetry of the matrices given by the mixed 
variational principles is not effectively utilized in the solution procedure 
using the IMSL routine EIGZF (ref. 36). If the symmetric nature of these 
matrices were effectively used, the computation time might be further reduced. 
Similarly, no advantage has been taken of the banded nature of the matrices in 
using the finite-difference methods. If this advantage were utilized in the 
computational scheme, some CPU time and considerable memory size might be 
saved. 

Finally, the mathematical implications of the various analytical methods 
are considered. 

In the two finite-difference procedures presented, the formulation of the 
eigenvalue problem does not require great mathematical skills or judgment. 
Anyone with a reasonable mathematical background can formulate and solve the 
problem, once the recursive relations for elimination of the fictitious 
stations are provided. These methods require no integrations and thus the 
necessity of using any numerical quadrature formulae is avoided. For the class 
of problems treated so far, no numerical instabilities were observed, meaning 
that the finite-difference procedures lead to numerically well-posed problems. 
Further, for untwisted blades, the coupled or uncoupled frequencies converge 
from below and close lower bound solutions are obtained. 

Considering the mixed variational principles, one starts with a mathemat- 
ically stationary functional (which possesses neither a maximum nor a minimum) 
rather than a minimum energy concept. Consequently, the natural frequencies 
do not specify any bound but oscillate closely around a converged exact value. 
The convergence rates and accuracies depend largely on the choice of shape 
functions. The final equations require evaluation of definite integrals using 
a numerical quadrature formula generally. Thus, some judgment and skill of 
the investigator are needed when the mixed variational principles are used. 
Furthermore, the matrices are usually ill-conditioned, and some form of scaling 
is necessary if one wishes to evaluate the frequencies by determinant search. 
While the mass and stiffness matrices are symmetric, they are not necessarily 
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nonsingular. One may thus have computational problems if inversions are 
needed. These problems could however be obviated by condensation or by rewrit- 
ing the equations so that the elements along the leading diagonal are nonzero, 
thus losing the symmetry condition. Once the mixed variational principle is 
formulated and suitable shape functions are assumed, the solution procedure is 
rather straightforward however. 

In conclusion, each class of methods has certain distinct advantages and 
one can choose a method depending on the constraints at hand, considering, the 
advantages and disadvantages. All the methods presented have capabilities of 
producing accuracies of practical interest. The second order finite-difference 
procedure is definitely superior to the classical approach while for the class 
of problems considered here, the Reissner method is superior to the Dean and 
Plass principle. 


CONCLUDING REMARKS 

The coupled bending-bending-torsion vibrations of pretwisted blades of 
asymmetric airfoil cross section are analyzed by solving the governing coupled 
equations of motion using first and second-order finite-difference methods. 
Alternative approaches of solving the coupled vibration problem are demonstra- 
ted by employing two mixed variational approaches, one proposed by Dean and 
Plass and the other by Reissner. The relative merits, convergence properties 
and accuracies of all the four methods are discussed. The effects of asymmetry 
and pretwist in coupling the component modes are studied. In the course of 
this study, the following conclusions have emerged: 

(1) The improved second order finite-difference procedure is far supe- 
rior to the classical first order procedure. One can use the present improved 
theory for attaining better accuracies while using smaller computation time 
and effort than that required by ordinary finite-difference procedure. 

(2) For untwisted blades performing coupled bending-bending-torsion 
vibrations, close lower bound solutions can be obtained by using improved 
finite-difference procedure. However pretwist disturbs the bound offered by 
finite-difference procedures. The mixed variational approaches do not indicate 
any specific bound. 

(3) When identical shape functions and beam kinematics are used, the 
Dean and Plass and the Reissner methods give identical results. However, the 
Dean and Plass method requires more computational time. 

(4) The Dean and Plass method appears to be more advantageous when gyro- 
scopic effects are included in the analysis. Further investigation in this 
direction appears to be necessary for utilizing the latent potential of this 
relatively new variational principle. 

(5) For general vibratory problems, the variational methods seem to be 
more appealing in terms of computational space and time. However, some care 
and judgment of the investigator are required when the mixed variational prin- 
ciples are employed. 
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The finite-difference method is relatively straight forward and simple. 
One can apply this approach in a routine manner once the recursive relations 
are given for elimination of fictitious stations. 

(6) The effects of asymmetry and pretwist appear to be difficult to 
understand at first sight. However, once the coupled frequencies are classi- 
fied as belonging to a basic mode category as described in this work, the cou- 
pling trends can easily be established. 
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APPENDIX A 


DEVELOPMENT OF THE DEAN AND PLASS VARIATIONAL FUNCTIONAL 

The Dean and Plass dynamic variational functional (ref. 24) can be 
written as 


- Iff icA - 


P.P. 

~2^ } dv " J RS 


(Al) 


where I R $ is the Reissner functional given by (ref. 22) 

I DC = JJJ' {x..e.. - U ( t . - ) } dv - Jjf B . U .dv - ff T . U .ds 

RS lj 1J 0 v 1 j' l i Vi 1 l 


(A2) 


SI 


Development of the Reissner functional is presented in reference 12 for the 
general case of asymmetric blade. Only the kinetic energy equivalent of the 
Dean and Plass dynamic functional will be developed in what follows. 

Neglecting the effects of shear deformation, the general displacement 
field for coupled bending-bending-torsional motion can be written as (ref. 12) 


U x = x - ye = x : - y : e 
U y = y + xe = y x + x x e 

U z - - *l x l - *1*1 + 0 c 9 ' 


(A2) 

(A3) 

(A4) 


For the computation of the kinetic energy, one can disregard the effects of 
longitudinal inertia effects in the total kinetic energy without great loss of 
accuracy in the final results. It will be assumed that the longitudinal 
inertia effects are negligible and thus the linear momentum densities about 
x and y directions take the following forms: 


\ - 0 fl x - P ( J 1 - *1 5 > 


(A5) 


V 


P u y = P (yi + 5J0) 


(A6) 


Equations (A5) and (A6) contain terms associated with linear velocities xj 

and yi together with the angular velocity e. One can explicitly define 
the linear and angular momenta by defining the linear momenta about x^x^, yjy^ 
centroidal directions and angular momentum due to elastic twist independently 
as follows: 



(A7) 
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P yl = pyl 


(A8) 


P 0 


plp9 

r~r t 


(A9 ) 


Thus 


d '> -| L •C|| P xl ■ p *l®> ( *1 " » 1 ®> + (P yl + p il®> (j'l * * 1 ®)] 

= f Q [ P xl A *l + P yl A >1 * p e A ® ] d2 


dA dz 


(A10) 


and 


rrr^^J 

~ 2 2 2 2 " 

P ,A + Pi A P^A^ 
xl .yl u. 9 

J JJ 2p «-J 0 

CL 

Q. 

CM 

a 

CM 

1 


dz 


(All) 


The Reissner functional, Irs> can be reduced to the following form for the 
coupled bending-bending-torsion case neglecting shear deformation, body forces, 
surface fractions and warping rigidity effects: 


*RS 




.. " ..." T . T e . "x’ylyl " 2M xVxlyl + "y'xlxl 
y 2E ('xlxl'ylyl - 'xV . 


dz 


(A12) 


Making use of equations (A10) to ( A12 ) in equation (Al), the Dean and Plass 
dynamic variational functional becomes 




/ 

0 


p2 a A P^A^ 

P xi A (x + r y 9) + P yl A (y + r x 0) + P 0 A0 - -|^ - -i, 


2p 2pT 


T ' . i . M xW - 2M xVxlyl * "y'xlxl 
+ M..y + M.x - T_0 + yrr- + 


' 0 W 7T 


2 E ^xlxlVyl " T xlyl^ 


dz (A13) 


Treating P x i, Pyi, P 0 , M x , My, T 0 , x, y and e as independent variables 
(independent of one another but dependent on z and t) and performing the varia- 
tional process (ref. 12) 


6 ( 1) In = 0 


(A14) 
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one obtains the following equations: 


xl 


= p 


(x + re) = P x 1 


\ 


p y i = p (y + r x e) = pyi 


P e = p ^xlxl + ^lyl^ ® /A = pI p ® /A 

y " = _ (M x ! ylyl ~ Vxlyl^ 

E ^xlx^ylyl " ! xlyl^ 


X as - 


<¥xlxl - Vxlyl* 

E ^xlx^ylyl " ! xlyl^ 


e =T 0 /2C 


J 


(A15) 


- It < AP xl> + 77 <V - 0 

a Z 


- It (AP - > 


yi' 7? l "x 

J 3Z 


(MJ = 0 


(A16) 


- (AP ,r + AP ,r + AP ) + — (T ) = 0 
at v xl y yl x e az e' 


— (M^) = 0 at z = 0, L or x prescribed 


— (M v ) = 0 at z = 0, L or y prescribed 

d Z X 


(A17) 


T = 0 at z = 0, L or a prescribed I 

/ 

Equations (A15) relate the linear and angular momenta to respective linear or 
angular velocities, curvatures to bending moments and rate of twist to twisting 
moment and torsional rigidity. These relations can alternatively be looked 
upon as stress-strain relations. Equations (A16) represent the conditions of 
equilibrium or motion. The first two of equations (A16) relate the time rate 
changes of linear momenta per unit length to corresponding rates of loading 
(forces per unit length). The last of equations (A16) equates the time rate of 
change of angular momentum per unit length to the twisting moment (torque) per 
unit length. Equations (A17) are the boundary conditions associated with the 
present formulation. Assuming harmonic solutions for eliminating the time 
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dependence, equations (A16) are rewritten in terms of x, y and e as 
follows: 


“7 j EI xlxl “^T + EI xlyl prj = pAP n y + pAP n r x e 


^-(ei &+EI ^2- 

dz 2 fyiyi dz 2 EI *iyi dz 2 


| . pAP 2 x 


+ ' AP nV 


dz (<#) - - ‘“'Vn* - pAr y P n x ‘ pI cf P n‘ 


(A18) 


Equations (A18) agree with those derived earlier by Carnegie (ref. 10) for the 
specialized case considered here. 
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APPENDIX B 


SHAPE FUNCTIONS, GEOMETRIC PROPERTIES, AND COEFFICIENTS 
OF DIFFERENTIAL EQUATIONS 

Shape functions: 

f i (n) = n 1 (Rl ti - R 2 n 2 + R 3 n 3 ) 

9i = R 4 ( 1— n ) 1 
hi = n i (1 - R 5 n) 

G(t) = sin Pt 

H(t) = cos Pt 

Rj. = (i+2)(i+3)/6; R 2 = i(i+3)/3; R 3 = i(i+l)/6; R 4 = l/(i+l); R 5 = i/(i+l) 

Geometric properties: 
c x = r xo Yn — c y gSi n yn 

r y = r y0 cos yn - r xo sin Yn 

Ixlxl = lyy sin 2 Yn + I xx cos 2 yti 
lylyl = lyy cos 2 Yn + Ixx sin 2 Yn 



Coefficients of differential equations: 

The coefficients of the differential equations given by equations (10 to 12) 
are given below. 

a l = x yy s i n ^Yn + Ixx cos^Yn 

bi = 2y sin 2 yh (lyy - I X x) 

Cl = 2y 2 cos 2yn (lyy - Ixx) 

di = ( lyy — Ixx) sin 2 Yn/ 2 

ei = 2 y cos 2 yp (lyy - Ixx) 

fl = -2 y 2 sin 2 yt> (lyy - I xx ) 

a 2 = Uyy - Ixx) sin 2 Yn/2 

b 2 = 2 y COS 2yT 1 (lyy - I xx ) 
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C2 = -2^r 2 sin 2yn (Iyy - Ixx) 
d 2 = Iyy cos^ Y n + Ixx s '* n ^ Yn 
e 2 = -2y sin 2yn (Iyy - Ixx) 
f 2 = -2y 2 cos 2yn (Iyy - Ixx) 

9 ! = P AL 4 /E 

hi = pAL 4 (r x0 cos yn - r y0 sin yn)/E 
92 = 91 

h2 = pAL 4 (r y0 cos yn + r x0 sin yn)/E 

313 = 1 ; d3 = p [I xx + Iyy + A (r £ 0 + r^ 0 )]/C 

b3 = pA (r x0 cos yn - r y0 sin yn)/C; n = i h 

C3 = pA (r y0 cos yn + r x0 sin yn)/C; h = L/n 
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TABLE I. - NUMERICAL DATA FOR AIRFOIL 
BLADES (ref. 30) 


L = 152.40 mm (6.0 in. ) 

A = 589.68 sq mm (0.914 sq in.) 

I xx = 34.963 mm 4 (0.000084 in, 4 ) 

I yy = 2792.9 mm 4 (0.00671 in. 4 ) 

E = 213.73747 GPa (31xl0 6 psi) _ 
p = 20.3444 kg/rrr ( 0.284 lbf/in.) 
r XQ = 0.19304 mm (0.0076 in.) 
r yo = 1.1938 mm (0.047 in.) 


TABLE II. - TORSIONAL RIGIDITY OF PRETWISTED 


AIRFOIL BLADE (ref. 29) 


Pretwist angle, 

r, 

deg 

Torsional rigidity, C 

Nm^/rad 

lb in. 2 /rad 

0 

9.1265 

3180.1684 

15 

9.1759 

3197.4 

30 

9.3269 

3250.0 

45 

9.6282 

3355.0 

60 

10.0440 

3500.0 

75 

10.6470 

3710.0 

89 

11.5600 

4028.0 


TABLE III. - CONVERGENCE PATTERN ANO COMPARISON OF COUPLEO BENDING-BENDING-TORSION FREQUENCIES: FIRST- AND 


SECOND OROER FINITE-OIFFERENCE METHODS (Hz) 


■ 

Mode 

number 

First order finite-difference method 

Second order finite-difference method 

n . 10 

20 

30 

40 

50 


10 

20 

30 

Extrapolated 


M 

95.929 

96.534 

96.648 

96.687 


96.716 

96.706 

96.734 

96.737 

96.738 



581.449 


603.358 

604.613 

605.196 



605.911 

606.141 

606.204 

IWI 

n 

836.642 

841.330 

842.204 

842.510 

842.652 

842.729 

842.658 

842.874 

842.895 

842.901 



1074.560 

1076.166 

1076.468 

1076.574 


EHSEil 


E5ZSZI 

1076.709 

1076.710 


■Jl 

1557.277 

1659.921 

1680.532 

1687.877 


1693.168 

1672.018 

1694.623 

1696.647 

1697.190 


6 

2875.489 

3087.379 

3091.582 

3093.051 


EZHEJ 

3093.880 

3094.862 

3094.923 

3094.937 

■ 

7 

3064.882 

3201.482 

3269.55 2 

3294.057 

3305.523 

3311.786 

3207.516 

3313.419 

3322.705 

3325.192 


np 

103.644 

98.249 

98.552 

97.324 

97.221 

97.165 

101.665 

98.171 

97.554 

97.378 



384.472 

482.100 

501.378 

508.273 

511.493 

513.251 

464.308 


513.490 

515.149 



1032.234 

988.711 

980.034 

976.911 

975.446 

974.645 

1011.948 

981.844 

976.653 

975.145 



1187.743 

1109.888 

1104.447 

1102.855 

1102.162 

1101.796 

1134.652 


1102.703 

1101.983 



1187.743 

1548.318 

1607.632 

1628.241 

1637.768 

1642.942 

1490.641 

1632.749 

1647.241 

1651.192 


aTfO 

2622.674 

3089.437 

3129.428 

3133.785 

3135.325 

3136.082 


3139.478 

3138.599 

3138.141 


u 

3112.428 

3151.108 

3213.081 

3245.109 






3290.033 



























TABLE IV - COVERGENCE PATTERN OF COUPLED BENDING-BENDING- 
TORSION FREQUENCIES (Hz) 


(a) Reissner method and Dean and Plass method using 
identical shape functions, y = 0° 


Mode 

number 


Number of terms in 

assumed solution, k 


i 

2 

3 

4 

5 

6 

n 

94.042 

96.664 

96.738 

96.739 

96.739 

96.739 


822.909 

573.170 

599.444 

606.259 

606.219 

606.235 

El 1 

1079.098 

842.363 

842.903 

842.904 

842.904 

842.904 



1076.708 

1076.708 

1076.711 

1076.711 

1076.711 



3183.905 

1596.621 

1629.061 

1697.967 

1696.414 

B 


5216.526 

3107.792 

5293.266 

3095.274 

3210.000 

3051.492 

3095.062 

3094.935 

3328.162 


(b) Reissner method and Dean and Plass method, Y = -30‘ 






















TABLE V. - NATURAL FREQUENCIES OF DOUBLY SYMMETRIC PRETWISTED BLADES (Hz): 


SECOND ORDER FINITE-DIFFERENCE METHOD n = 30, r xo = r y0 = 0 


Mode 

number, 

classification 

Uncoupled 
bending and 
torsion modes 
r xo = r yo = y = 0 

Coupled bending-bending modes and uncoupled 
torsion modes 

y = 30° 

II 

c_n 

0 

Y = 60° 

II 

00 

uo 

o 

1 

(FI) 

96.738 

97.558 

98.417 

99.319 

100.255 

2 

(F2) 

606.158 

513.908 

454.448 

403.476 

330.908 

3 

(SI) 

864.603 

1014.142 

1136.247 

1250.120 

1282.718 

4 

(Tl) 

1051.080 

a 1062. 557 

a l079.585 

a 1102.668 

a 1182.920 

5 

(F3) 

1696.742 

1648.125 

1600.250 

1562.417 

1674.196 

6 

(T2) 

3153.229 

a 3187 . 661 

a 3238 .745 

a 3307. 992 

a 3548.748 

7 

( F4) 

3323.060 

3281.225 

3234.039 

3176.773 

3071.047 

8 

(T3) 

5255.263 

a 5312 . 648 

a 5397.786 

a 5513 . 196 

a 5914.446 

.9 

(S2) 

5417.605 

5242.015 

5149.207 

5047.868 

4815.582 

10 

(F5) 

5488.150 

5815.830 

6086.129 

6417.764 

7315.336 


Uncoupled torsional frequency of pretwisted blade. 


TABLE VI. - COUPLED FREQUENCIES (Hz.) OF PRETWISTED BLADES WITH 


ASYMMETRY ABOUT ONE PLANE: r x0 = 0.193 mm, r yo = 0 


Mode 

number 

Uncoupled 
frequency 
r xo = r yo = y = 0 

Coupled 

bending- 

torsion 

frequency, 

Y = 0 

Coupled bending-bending-torsion 
frequency 

y = 30° 

-< 

ii 

0 

Y = 60° 

Y = 89° 

wm 

96.738 (FI) 

96.737 

97.558 

98.417 

99.319 

100.255 


606.158 (F2) 

606.142 

513.891 

454.432 

403.461 

330.898 


864.603 (SI) 

a 864.603 

1013.977 

1136.514 

1250.296 

1283.036 

IBM 

1051.080 (Tl) 

1051.095 

1062.743 

1079.348 

1102.544 

1182.623 


1696.742 (F3) 

1696.647 

1648.040 

1600.157 

1562.299 

1674.142 

Hi 

3153.229 (T2) 

3153.270 

3187.303 

3245.592 

3309.104 

3549.324 

n 

3323.060 (F4) 

3322.728 

3281.327 

3226.973 

3175.471 

3070.251 


Uncoupled flexural mode in chordwise direction 























TABLE VII. - COUPLED FREQUENCIES (Hz) OF PRETWISTED BLADES WITH 
ASYMMETRY ABOUT ONE PLANE: r xo = 0, r y0 = 1.1938 mm 


Mode 

number 

Uncoupled 

frequency...;..-.;. 

r xo = r yo = y = 0 

Coupled 

•.bending- 

torsion 

frequency, 

Y = 0 

Coupled bending-bending-torsion 
frequency 

o 

O 

oo 

II 

Y = 45° 

0 

O 

<£> 

II 

-< 

ii 

CO 

VO 

0 

■■1 

96.738 (FI) 

a 96.738 

97.556 

98.412 

99.311 

100.246 


606.158 (F2) 

a 606.158 

513.342 

454.127 

403.354 

330.907 

El . 

864.603 (SI) 

842.893 

974.417 

1170.423 

1264.208 

1289.706 


1051.080 (Tl) 

1076.696 

1105.563 

1047.110 

1088.446 

1174.069 


1696.742 (F3) 

a 1696.742 

1647.418 

1599.494 

1561.908 

1671.711 

! — 

3153.229 (T2) 

3094.865 

3131.506 

3270.061 

3304.629 

3542.868 

mm 

3323.060 (F4) 

a 3323 .060 

3290.388 

3166.320 

3154.937 

3069.589 


Uncoupled flexural mode in flapwise direction. 


TABLE VIII. - COUPLED BENDING-BENDING-TORSION FREQUENCIES OF ASYMMETRIC 
AIRFOIL BLADE (Hz): r x0 = 0.193 mm, r y0 = 1.1938 mm 


















TABLE IX. - COMPARISON OF CPU TIMES (milliseconds) REQUIRED BY DIFFERENT 
METHODS FOR PRETWISTED AIRFOIL BLADE CASE USING EIGZF ROUTINE: 


y = -30 , r x0 = 0.193 mm, r v0 = 1.1938 mm 


yo 


Number 
of beam 
segments 
n 

Finite-difference method 

Number of 
terms in 
assumed 
solutions 
K 

Mixed variational principle 

I order 

II order 

Reissner 

Dean and PI ass 

10 

831 

856 

1 

59 

84 

20 

5 550 

5390 

2 

215 

348 

30 

17 736 

17014 

3 

490 

847 

40 

41 484 


4 

924 

1639 

50 

76 966 


5 

1494 

2881 

60 

225 838 

, 

6 

i 

2193 

4410 
















NATURAL FREQUENCY, Hz 



PRESENT RESULTS (SECOND- 
ORDER FINITE DIFFERENCE 
METHOD, REISSNER METHOD 
OR DEAN AND PLASS METHOD) 

CARNEGIE, etal. 130) 


(a) First, second, third, and fourth coupled mode frequencies. 



(b) Fifth, sixth, and seventh coupled mode frequencies. 

Figure 2. - Comparison of theoretical and experimental coupled 
bending-bending-torsion frequencies. 



RELATIVE AMPLITUDE 


RELATIVE 



(a) First coupled bending-bending mode lc) First uncoupled torsion mode 

(10a 255 Hz). (1182. 920 Hz). 



AXIAL FRACTIONAL LENGTH 


(b) Second coupled bending-bending mode (d) Third coupled bending-bending mode 

(33a 908 Hz). (1282.718 Hz). 

Figure 3. - Coupled bending-bending, and uncoupled torsion mode shapes of a pretwisted blade of 
doubly symmetric cross section. r XQ - r y0 • 0, y ■ 89°. 




(e) Fourth coupled bending-bending mode 
(1674196 Hz). 



(3071.047 Hz). 



(g) Second uncoupled torsion mode 
(3548.748 Hz). 


Figure 3. - Concluded. 


RELATIVE AMPLITUDE 




(a) First coupled bending-bending-torsion (c) Third coupled bending-bending-torsion 
mode (100.246 Hz). mode (1171910 Hz). 



AXIAL FRACTIONAL LENGTH 


(b) Second coupled bending-bending-torsion (d) Fourth coupled bending-bending-torsion 
mode (330. 896 Hz). mode (1291 688 Hz). 

Figure 4. - Coupled bending-bending-torsion mode shapes of a pretwisted asymmetric airfoil 
blade. r x0 - 0. 19304 mm (0. 0076 in), r y0 ■ 1. 1938 mm (0. 047 in), y • 89°. 



(e) Fifth coupled bending-bending-torsion 
mode (1672. 784 Hz). 



(0 Sixth coupled bending-bending-torsion 
mode (3067. 755 Hz). 



(g) Seventh coupled bending-bending-torsion 
mode (3539. 662 Hz). 


Figure 4. - Concluded. 
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